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ABSTRACT 

Many  algorithms  for  polynomial  least  squares  approximation  of  real-valued  function  on  a 
real  interval  determine  polynomials  that  are  orthogonal  with  respect  to  a  suitable  inner  prod- 
uct defined  on  this  interval.  Analogously,  it  is  convenient  to  compute  Szego  polynomials,  i.e., 
polynomials  that  are  orthogonal  with  respect  to  an  inner  product  on  the  unit  circle,  when  ap- 
proximating a  complex-valued  function  on  the  unit  circle  in  the  least  squares  sense.  It  may  also 
be  appropriate  to  determine  Szego  polynomials  in  algorithms  for  least  squares  approximation 
of  real-valued  periodic  functions  by  trigonometric  polynomials.  This  paper  is  concerned  with 
Szego  polynomials  that  are  defined  by  a  discrete  inner  product  on  the  unit  circle.  We  present 
a  scheme  for  downdating  the  Szego  polynomials  and  given  least  squares  approximant  when  a 
node  is  deleted  from  the  inner  product.  Our  scheme  uses  the  QR  algorithm  for  unitary  upper 
Ilessenberg  matrices.  We  describe  a  data-fitting  application  that  illustrates  how  our  scheme  can 
be  combined  with  the  fast  Fourier  transform  algorithm  when  the  given  nodes  are  not  equidistant. 
Application  to  sliding  windows  is  discussed  also. 


Downdating  of  Szego  Polynomials 
and  Data  Fitting  Applications 

Cl.S.  Amni.ii"      W.B.  Giagg*       L.  Reichcl: 


Abstract 

Many  algorithms  for  polynomial  least  squares  approximation  of  a  real-valued  function 
on  a  real  interval  determine  polynomials  that  are  orthogonal  with  respect  to  a  suitable  inner 
product  defined  on  tins  interval.  Analogously,  it  is  convenient  to  compute  Szego  polynomials, 
i.e..  polynomials  that  are  orthogonal  with  respect  to  an  inner  product  on  the  unit  circle, 
when  approximating  a  complex-valued  function  on  the  unit  circle  in  the  least  squares  sense. 
It  may  also  be  appropriate  to  determine  Szego  polynomials  in  algorithms  for  least  squares 
approximation  of  real-valued  periodic  functions  by  trigonometric  polynomials.  This  paper 
is  concerned  with  Szego  polynomials  that  are  defined  by  a  discrete  inner  product  on  the 
unit  circle.  \\<  present  a  scheme  tor  downdating  the  Szego  polynomials  and  given  least 
squares  approxiinaut  when  a  node  is  deleted  from  the  inner  product.  Our  scheme  uses  the 
Oil  algorithm  lor  unitary  upper  llessenbei  ^  matrices.  We  describe  a  data-fitting  application 
that  illustrates  how  our  scheme  can  be  combined  with  the  fast  Fourier  transform  algorithm 
when  the  given  nodes  are  not  equidistant.  Application  to  -luling  windows  is  discussed  also. 


1      Introduction 

Let  {-/c}"l=l  he1  a  set  of  distinct  nodes  on  the  unit  circle  and  h*t  {it'?}™^  be  n  set  of  positive 
weights.  Introduce  for  complex- valued  functions  g  and  h  defined  .ii  i  he  nodes  the  discrete  inner 
product  on  the  unit  circle 

(uJ')m  :=  ^2'A=kMsk)wl,  (1.1) 

A  =  l 

where  the  bar  denotes  complex  conjugation.  Polynomials  that  are  orthogonal  with  respect  to 
an  inner  product  on  the  unit  rirrlo  are  known  as  S:t  go  polynomials.  Let  {pj}''!"1  denote  the 
family  of  ortlionornwl  Szego  polynomials  with  respect  to  the  inner  product  (1.1),  where  dj  is  of 
degree  j  and  has  a  positive  leading  coefficient.  The  polynomials  <?>_,  satisfy  the  Szego  recurrence 
rr  la  lions 

Po(-)     =     Oo(~)  =  1/ffo, 
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rrJ  +  loJ  +  n:)     =     =oj(z)  +  ';j+iQj{z),    0  <  j  <  m-1,  (1.2) 

aJ+idj+i(z)     =     s"fj+i<i>j{z)  +  dj(z), 

where  the  recurrence  coefficients  7J  +  1  t  C  and  a/  +  1  >  0  are  determined  by 

7;+i     =     -d.coj,,, /<\.  (1.3) 

/  ,.,\l/2 

^,  +  i     =     (1  -  iTj+iTJ       •  0  <  j  <  m, 

See,  for  example,  Grenander  and  Szego  [L  Chapter  2].  It  can  be  shown  by  induction  that  the 
auxiliary  polynomials  o;  in  (1.2)  satisfy 

Oj(  =  )  =  =Jd>j(l/:).      ()<j<m, 

where  d>j(z)  denotes  the  polynomial  obtained  by  conjugating  the  coefficients  of  6j{z)  in  the 
power  basis.  Since  the  measure  on  the  unit  circle  that  defines  (1.1)  has  precisely  m  points  of 
increase,  we  have  |-.;|  <  1  for  1  <  j  <  in  and  |7m|  =  1.  The  coefficients  "j  are  known  as  Schur 
parameters,  and  we  refer  to  the  a1  as  the  associated  complementary  parameters .  Although  the 
complementary  parameters  are  mathematically  redundant,  we  retain  them  during  computations 
in  order  to  avoid  numerical  instability.  Note  from  (1.3)  that  a^  is  the  total  weight  of  the  measure 
that  defines  (1.1),  and  that  6~  is  the  leading  coefficient  of  <t>j  for  0  <  j  <  m. 
For  later  reference  we  also  define  the  polynomial 

Om{z)  :=  :om-,(:)  +  1m4>m-\.{*)  -  (1-4) 

Then 

<U*)  =  C-i  H(z-zk),  (1.5) 

Jk=l 

In  particular,  (0m,cV,)  =  0  for  0  <  j  <  m. 

Example  1.1.  Let  Zk  :=  exp(2/T7(A:  —  1 )/???)  and  w\  :=  1  for  1  <  k  <  m.  where  i  :=  \J  —  1. 
Then  oq  =  m1'2,  o3  =  1  and  yj  =  0  for  1  <  j  <  m,  and  -ym  =  I.  Thus.  d>j(z)  =  ?n-1//2cJ  for 
0  <  j  <  m,  and  Om{  =  )  =  mTxl2{zm  -  1).  □ 

Introduce  the  discrete  norm 

\\g\\m  :=  (fl.fl)!/2, 


and  let  II„_i  denote  the  set  of  all  polynomials  ol  degree  at  most  n  -  1.  Let  /  be  a  given 
complex-valued  function  defined  at  the  nodes  r,  .  and  consider  the  problem  of  computing  the 
polynomial  pn-\  t  Hn-i-  for  some  n  <  in.  that  satisfies 

11/  -  Pn-\\\m  =     iiijn     \\f  -  p\\m  •  (1.6) 

The  solution  pn-i  of  the  minimization  problem  (  l.(i)  can  be  expressed  conveniently  in  terms  of 
the  Szego  polynomials  o:  as  follows.  Introduce  the  vector 

f  =  ["i/Ci  ).«-,./  •(.:.  I ir,nf[sm)]T  .  (1.7) 

and  the  in  x  u  matrix  O  — 

<lkj  :=  Oj-X[zk)trk,         _       <  m.     1   <j  <  n.  (1.8) 

where  t/u.  :  =  \/wr.    Note  that  the  matrix  O  has  oil  honormal  columns,  i.e..  O'O  —   I.  where 

V       k  •  -       - 

C)":=QI.   Let  the  vector  a  =  [n-i.o-j m,,]'    be  delined  by 

a:=C;wf.  (1.9) 

Then  the  polynomial  />u_i  can  be  written  as 

/,„_,    =   yit,o ,_,  .  (1.10) 

=  I 

where  the  Fourier  coefficients  n7  arc  independent  of  a. 

It  is  the  purpose  of  this  paper  to  present  an  algorithm  for  downdating  the  recurrence  coef- 
ficients 7j  and  tr,.  as  well  as  the  Fourier  coefficients  n  .  when  an  node-weight  pair  is  removed 
from  the  inner  product  (1.1).  Our  algorithm  is  based  on  i  he  observation  that  the  columns  of  the 
matrix  0  arc  eigenvectors  of  a  unitary  IfcssenlxM",  matrix  defined  by  the  recurrence  relations 
( 1.2)— ( 1.3).  This  makes  it  possible  to  downdate  ilie  coefficients  ~(j,  a3  and  ctj  by  applying  the 
QR  algorithm  for  unitary  llesscuberg  matrices,  presented  in  [o],  with  the  node  to  be  removed 
as  shift.  Details  are  described  in  Section  2. 

We  remark  that  the  problem  of  updating  i  he  coefficients  ~fj,  (jj  and  aj  when  a  node-weight 
pair  {zm+i,  ".'4  +  1 }  's  added  to  the  inner  product  ( 1.1 )  is  discussed  in  [8].  The  updating  problem 
can  be  solved  by  using  an  inverse  QR  algorithm  for  unitary  Hessenberg  matrices:  see  [8,  1]. 

Assume  that  we  wish  to  determine  i  he  polynomial  pn-\.  given  by  (1.0 )— ( 1.10)  with  n  <  m 
when  the  set  of  m  nodes  in  the  inner  product   (1.1)  is  a  subset  of  the  set  of  .V  equidistant 

:$ 


nodes  {exp(2;r/j/.Y  j}'^1-  with  «  :=  A  -  in  >  0  small.  The  weights  wj:  are  all  assumed  to  be 
unity.  Then  it  may  be  attractive  to  compute  pn-\  by  firs)  computing  the  polynomial  interpolant 
/>v-i  €  n.v-i  on 'the  set  ol"  .V  equidistant  nodes  by  lining  the  fast  Fourier  transform  (FFT) 
algorithm,  and  then  determining  pm~\  from  /'.y-i  by  applying  our  downdating  scheme.  The 
Fourier  coefficients  of  the  polynomial  approximant  pn-\  are  then  equal  to  the  first  n  Fourier 
coefficients  of  pm-\-  This  approach  can  similarly  be  applied  to  trigonometric  polynomials. 
Details  are  described  in  Section  3.  Computed  examples  are  presented  in  Section  4  and  Section  5 
contains  a  summary. 

Updating  and  downdating  of  polynomials  approximants  /;n-i  when  all  the  nodes  Zk  are 
real  has  received  considerable  attention  in  the  literature:  see  Scott  and  Scott  [9]  and  references 
there.  A  collection  of  algorithms  for  updating  and  downdating  based  on  orthogonal  polynomials 
is  presented  by  Flhay  et  al.  [3].  Algorithms  for  updating  are  also  considered  in  [G,  "]. 

2      Downdating  of  Szego  polynomials 

The  connection  between  the  Szego  polynomials  determined  by  (1.1)  and  a  unitary  llessen- 
berg  matrix  can  be  seen  as  follows.  Using  the  basis  of  Szego  polynomials,  we  can  write 

J2>lijOi-[(:)  =  :Oj-\(;).    1  <  j  <  m.  (2.1) 

where  J/j+i.j  —  a  i  '01'  1  ^  J  <  "/.and  ij,„  +  lm  =  I.  l.el  ///. ,  :=  0  for  3  <  j+2  <  k  <  m,  and  define 
the  upper  Ilesscnberg  matrix  //  =  ['M/]"V=i  •  Also  (hTme  the  unitary  matrix  U  =  [/U-jlTfc- 1  ^y 

fikj  ■=  Oj-i(^.)'a  •  (2.2) 

Substitution  of  z  —  z^,  1  <  k  <  in,  into  (2.1)  yields  the  equation 

UTUT  =  UT\,  (2.3) 

or.  equivalently, 

//  =  ("Ac/.  (2.4) 

where  A  :=    diag[ci.~j zm].   Thus.  If  is  a  unitary  upper  Ilesscnberg  matrix  with  positive 

subdiagonal  elements  that  is  uniquely  determined  by  the  inner  product  (1.1).  Also  note  that 
the  first  n  columns  of  U  make  up  the  matrix  Q  define<l  by  ( 1  .S).  Algorithms  for  the  solution  of 
the  least-squares  problem  (  1.(5)  can  therefore  lie  viewed  in  terms  of  the  spectral  decomposition 

(2,1). 


It  is  fairly  straightforward  to  show,  by  using  the  recurrence  relations  ( 1.2)— (1.3),  that  // 
can  be  written  as  a  product  of  in  elementary  unitary  transformations  that  are  defined  by  the 
recursion  coefficients  -• ,  and  a,  for  the  or  We  have 


H   -  t'll"  1  )(.'l{~,2)  •  •  •^'m-l(7m-l)G,,m(7m)i 

where  the  Gj{jj).  1  <  j  <  m.  are  m  x  m  Givens  matrices 


(2.5) 


^(7i):  = 


.?        'j 


fm-j-] 


Hid  6"m(7„, )  is  the  diagonal  matrix 


0      -     ;  :=    diag[l.l l.-7m]- 

We  refer  to  the  representation  (2.5)  as  the  Schur  parametric  form  of  //. 

The  development  of  efficient  algorithms  for  eigenproblcms  for  unitary  Hessenberg  matrices  is 
facilitated  by  the  fact  that  every  unitary  upper  Hessenberg  matrix  with  nonncgative  subdiagonal 
elements  has  a  unique  Schur  parameterization,  lor  example,  when  the  implicitly  shifted  QR 
algorithm  is  applied  io  find  the  spectral  decomposition  of  a  unitary  Hessenberg  matrix,  a 
sequence  of  intermediate  unitary  Hessenberg  matrices  is  generated  that  converges  to  a  diagonal 
matrix.  In  [5]  the  QR  algorithm  for  unitary  Hessenberg  matrices  is  formulated  in  terms  of  the 
Schur  parameters  of  the  intermediate  matrices.  This  results  in  an  implementation  that  requires 
only  O(m)  arithmetic  operations  per  iteration  on  a  matrix  of  order  m. 

Assume  for  the  moment  that  the  matrix  //  and  scaling  factor  <7q.  are  given,  but  that  the 
nodes  z^  and  weights  w%  of  the  inner  product  (1.1)  are  not  explicitly  known.  Then  it  follows 
from  (2.4)  and  (2.2)  that  the  nodes  r/,  are  the  eigenvalues  of  //,  while  the  normalized  weight 
Wkl '<To  is  equal  to  the  first  component  of  the  normalized  eigenvector  corresponding  with  z^.  (We 
assume  that  each  eigenvector  is  scaled  to  have  unit  Kuclidean  length  and  a  nonncgative  first 
component.)  The  unitary  Hessenberg  QR  algorithm  can  be  used  to  compute  the  matrix  A  and 
vector  Ue\.  without  computing  the  rest  of  V.  with  not  more  than  0{m)  arithmetic  operations 
per  iteration.  Throughout  this  paper  e,  denotes  the  /'''  axis  vector  in  C"\  The  QR  algorithm 
determines  one  node-weight  pair  at  a  lime,  and  lor  each  pair  computed  the  order  of  the  unitary 


upper  Hessenberg  matrix  is  reduced  by  one.  so  that  the  reduced  Hessenberg  matrix  corresponds 
with  the  node- weight  pairs  that  have  not  yet  been  determined. 

We  remark  that  in  the  case  that  the  nodes  of  rite  discrete  inner  product  are  real,  then  the 
analogue  of  the  matrix  //  is  a  real  symmetric  tridiagonal  matrix  T  with  positive  subdiagonal 
elements.  This  matrix  T  contains  recurrence  coefficients  for  orthonormal  polynomials  that 
satisfy  a  three-term  recurrence  relation.  Golub  and  Welsch  [4]  proposed  the  use  of  the  QR 
algorithm  for  symmetric  tridiagonal  matrices  for  the  computation  of  the  node-weight  pairs 
associated  with  T.  This  algorithm  also  determines  one  node-weight  pair  at  a  time,  and  reduces 
the  order  of  T  when  such  a  pair  has  been  found. 

Conversely,  the  construction  of  //   from  the  inner  product  (1.1)  can   be  regarded  as  an 

inverse  eigenvalue  problem.    In  particular,  given  the  matrix  A  =  diag[^i,co rm]  and  the 

vector  qo  :=  <7q  [wi,wo,. . . ,  wm]  ,  we  can  perform  a  sequence  of  elementary  unitary  similarity 
transformations  whose  composition  results  in  an  m  <  in  unitary  matrix  U  such  that  the  matrix 


on  the  right-hand 


1      0 
0      I 


T  1 


*      qu 
qu       A 


1      07' 

o     r 


T 


2.6) 


is  of  upper  Hessenberg  form  with  positive  subdiagonal  elements.  (The  •  represents  an  arbitrary 
scalar  that  remains  unchanged.)  In  other  words.  (  "qo  =  <To&i,  and  U"\U  is  r.  unitary  upper 
Hessenberg  matrix  //.    Consequently,  //  is  the  matrix  corresponding  with  the  inner  product 

(1.1). 

The  transformation  of  A  to  //  can  be  pi-i  formed  using  0(m2)  arithmetic  operations  with 
the  inverse  unitary  QR  (IL'QR)  algorithm  described  in  [l].  The  1UQR  algorithm  is  an  updating 
procedure  because  it  incorporates  node-weight  pairs  one  at  a  time.  After  the  j1"  stage  of  the 
algorithm  has  been  executed,  the  j  x  j  unitary  Hessenberg  matrix  corresponding  with  the  inner 
product  determined  by  the  first  j  nodes  and  weights  has  been  obtained. 

In  [8]  the  approximation  problem  (1.0)  is  solved  using  the  IUQR  algorithm  to  construct 
the  Schur  parameters  of  the  unitary  Hessenberg  matrix  //.  The  elementary  unitary  matrices 
are  accumulated  against  the  vector  f  during  the  algorithm  to  obtain  the  vector  uf  Fourier 
coclficients  a  =  /.'"f  without  explicitly  forming  U.  In  this  way  we  obtain  the  interpolating 
polynomial  pm~i  that  is  the  solution  of  (1.0)  with  n  =  m  in  0(m2)  arithmetic  operations.  Of 
course,  the  partial  sums  of  the  Fourier  expansion  of  the  interpolating  polynomial  yields  the 
solution  (1.10)  of  (1.0)  for  each  n  <  m.  Moreover,  if  one  is  interested  in  computing  only  />n-i 


for  n  <  in.  then  the  algorithm  can  be  curtailed  so  that  only  O(mn)  arithmetic  operations  are 
required  for  the  computation  of  the  parameter-  \':j}1J=l.  {<7j}"=0,  and  {o,}*=1.  See  [8]  for 
details. 

Now  assume  that  we  have  solved  the  least-squares  problem  (1.6)  with  ;/  =  in  by  the  ntethod 
described  in  (Si.  so  that  sets  oi  Schui  parameters  {",1  =i  a"d  of  complementary  parameters 
{(Tj  }"L0  corresponding  with  the  inner  product  I  1 . 1  ).  as  well  as  sets  of  Fourier  coefficients  {o^}"!.! 
of  the  interpolating  polynomial  pm-\  are  explicitly  known.  In  the  downdating  problem,  we  seek 
l o  m ilve  the  corresponding  least-squares  problem  wit  h  one  term  deleted  from  (1.1).  In  particular, 
let  i  be  the  node  that,  is  to  be  deleted  from  the  inner  product  (1.1)  and  let  w  be  the  square-root 
of  the  associated  weight.  In  order  to  -implify  some  formulas  that  follow,  wo  assume  without 
loss  of  generality  that  i  =  ;,„  and  w  =  ir,„. 

Introduce  the  inner  product 


(j.li), 


V 


i,  i  :.  \li  I  r.  ).'/•" 


(2.7) 


and  discrete  norm 


\fj\ln-]   ■=  i<JJ.j)n~-\- 


We  seek  the  polynomial  \>m-i  €  lTm_2  such  that 

\\f  ~  f>m-:\[,n-\    =        mill  -    !>  \m-\ 


2.s: 


Denote  the  family  ol  orthonormal  Szego  polynomials  with  positive  leading  coefRcient  associated 
with  (2.7)  by  {Oj}''lSj.  Also  let  {"';}'_~'  and  \n   }'_"'  denote  the  sets  of  recurrence  coefficients 

for  the  Oj.  and   let    A   :=  diag[ri ~-     -\]  :il|d   q(,  :=  rr~   [ir{ M'm-l] '  •  where  cr{)  :=  [n^  — 

lb')1'2  is  the  total  weight  of  the  measure  that  defines  (2.7).   Then  from  the  above  discussion, 
there  is  a  unique  unitary  Itcssenbcrg  matrix  //  and  a  unique  unitary  matrix  U  such  that 

c"qo  =  ^oei 

and 

//  =  U"\0,  (2.9) 

where  e_,  denotes  the  7      axis  vector  in  C"  "   .    Moreover,  the  recurrence  coefficienis  jj  and 

cTj  arc  the  Schur  parameters  and  complementary  parameters.  resp<<  lively,  of  //.    The  optimal 

polynomial  is  then  given  by 

- 1 


whore  the  vector  of  Fourier  coefficients  a   =    ;/m-".> dm_i]r  is  given   by  a   :=   t'*f  and 

f  :=[t£»l/(Zl) lim-l/(*n,-i)]r. 

Our  scheme  for  downdating  is  based  on  the  observation  that  II  and  L  can  bo  computed 
from  //  and  U  by  applying  one  step  of  the  QR  algorithm  with  ''ultimate1'  shift  i.  together  with 
some  permutation  similarity  transformations,  to  determine  a  unitary  matrix  II'  such  that 
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Then 


and 
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■■  Ii"a  .  (2.11) 

Thus,  the  downdated  Fourier  coefficients  d;.  are  obtained  by  accumulating  each  elementary 
unitary  transformation  against  a.  An  efficient  implementation  is  obtained  by  using  the  QR 
algorithm  for  unitary  Hessenberg  matrices  described  in  [»]. 

We  now  assume  that  i  :=  -/  for  some  /,  1  <  /  <  m,  and  let  w  :=  w\.  One  step  of  the  QR 
algorithm  with  shift  i  applied  to  the  matrix  //  determines  a  unitary  upper  Hessenberg  matrix 
V .  such  that 


//'  :=  \"  II  \ 


II      0 


0 


because  i  is  an  eigenvalue  of  //.  It  is  easy  to  see.  however,  that  the  QR  algorithm  applied  to  // 
will  not  yield  the  required  //.  because  the  vector  rruei  will  not  be  transformed  as  required  by 
(2.10).  On  the  other  hand,  an  RQ  algorithm  for  //.  in  which  the  transforming  matrix  V  would 
be  a  product  of  Givens  matrices  in  the  reverse  order  of  (2.12)  below,  would  yield  the  desired 
similarity  transformation. 

Instead  of  modifying  the  unitary  Hessenberg  QR  algorithm  of  [5]  to  perform  an  RQ  iteration, 
we  apply  the  QR  algorithm  to  the  unitary  Hessenberg  matrix  II p  :=  .III  J,  where  J  = 
[e,n,em_i , . . . .  ei  ]  is  the  reversal  nuiirir  ol f  order  in.  It  is  easily  seen  that  if  //  is  given  by  (2.5), 
then 

//p  =  G'1(7i)r/>(72)--.r;m_1(7m.1)G,rri(77n), 

where  xfj  :=  Jm-jlm  for  1  <  j  <  m  and  -~,„  :=  -„.  The  application  of  the  QR  algorithm  on 
llp  is  equivalent  with  that  of  the  RQ  algorithm  on  // .  One  iteration  of  the  QR  algorithm  with 


shift  i  applied  to  II H  generates  a  unitary  upper  llcssenberg  matrix 


1 '  =  G , ( tS,  )6':( d2 )  •  •  ■  O'm-i ( (\n- 1  )6"_ (<5 
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such  that 


r-//pT/  = 


//'    o 
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Moreover.  Sm  can  he  taken  to  he  an  arbitrary  uniinodular  number  because  deflation  has  taken 


place  [ol.  Let  r,  :=  (1 


on!/ 


denote  the  complementary  parameter  to  6j  for  1  <  j  <  m. 


Observe  that  only  the  last  two  components  of  W0em  are  nonzero,  and  that  \Sm\  =  1  can 
bo  chosen  so  that  these  two  components  arc  given  by 
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Finally,  we  transform  the  right-hand  side  of  (2.1:5)  by  similarity  using 


,  where  J  is 


./      0 
0r      1 

the  revers. d  matrix  of  order  in  —  1,  and  transpose  the  result.   In  this  way  we  transform  11  to 
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where 


11'  =  JY 
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,  and 


0'       1 

by  the  uniqueness  of  the  reduction.  //     =  //.  ^orm_!  =  a0,  and  |c"m_i|ao  =  w.   The  vector  of 

Fourier  coefficients  a  is  then  determined  from  (2.11)  by  applying  W*  to  a  incrementally. 

Observe  that  our  downdating  procedure  requires  knowledge  of  the  node  to  be  deleted  but  not 
of  the  corresponding  weight.  We  can  therefore  compare  the  computed  value  of  ib  to  the  actual 
value  il'i  in  order  to  assess  the  accuracy  of  the  computations.  If  w  is  not  close  to  wi,  then  the 
downdatcd  polynomials  6:  and  pm-2  might  not  be  accurate.  Another  accuracy  check  is  provided 
by  the  computed  value  of  (>,„  in  the  algorithm.  This  quantity  is  the  m  '  diagonal  element  of 
the  upper  triangular  matrix  in  the  QR  factorization  ol  II1'  -  zl,  which  is  mathematically  zero 
when  z  is  an  eigenvalue  of  //.  If  the  computed  value  of  pm  is  not  "tiny",  then  the  computed 
downdatcd  polynomial  might  be  inaccurate.  We  therefore  consider  pm  and  ib  as  being  part  of 
the  output  of  the  algorithm. 


Algorithm  2.1  (Downdating  by  removing  the  node-weight  pair  {z.ib}  ) 
Input:   m.  ft;}*!,  fo}£0.  {Qj)?=l>  =  = 
Output:   ft;}^,  {^Y;=o-  {ojYjlV-P'n,  *: 

hA?J  ■=  TmlTm-ilSTi1;  E^il^Ti1  :=  k«-i]7=il; 
[Qil5Li:=[om-j+i];,l«i; 

60  ;=  £0  :=  l; 

for  j  :=  1,2 m  -  1  do 

/   2  -       2\i/2 

A'.;  :=  (^J  +  F«j-l  +  7j*j-l|  J        •' 

if  j  >  2  then  <7j_i  :=  Tj^ypji 

Tj  :=  ^-/pj-; 

<!>,  :=  (i^_!  +-rjfij-\)/pj; 

t\j  :=  -(^Oj  +  V1j  +  i; 

a,  +  1   :=  7-/.1,  +  'V-tj+l' 

rfj  :=  (^_i  +sfjz6j-l)/pj; 

Pm  :=   l-^77i-l  +  777i^m-l|/ 

<7m-i  :=  Tm_ipm;  w  —  I »Sn _  1 1 <t0 ,•  <t0  =  rm_iao; 


/m—  1    •  —    ,r7i  —  1 , 


fm-il/  <^?n-i  :—  0;  (parameter  correction] 


r  -  i77i  —  2  .       -  ri  im- 2   .       r  -  lm  — 2  .       r-  im- 2. 

[Tjjjsl     :=7rn-l[7m-j-lJj=l      ■        L^'Jj=l     •  =  F»>  -j  - 1  J7=  I    ■ 

[&j]?Jil--=l*m-j]TJix;  ° 

The  algorithm  overwrites  the  Fourier  coefficients  {'i,}'rL1  with  intermediate  quantities.  Al- 
gorithm  2.1  requires  0{m)  arithmetic  operations  (  +  .—.*,/)  and  the  evaluation  of  m  —  1  square 
roots. 

We  have  already  noted  that  if  the  nodes  Zk  are  real,  tlien  the  analogue  of  the  matrix  // 
is  a  real  symmetric  tridiagonal  matrix  T  with  positive  subdiagonal  elements.  Similarly  with 
Algorithm  2.1,  downdating  of  the  orthonormal  polynomials  associated  with  T,  as  well  as  of 
Fourier  coefficients,  can  be  carried  out  by  algorithms  based  on  the  QR  algorithm  for  symmetric 
tridiagonal  matrices.    This  observation  may  be  new. 

In  certain  data-fitting  applications  it  may  be  desirable  to  update  the  polynomial  pm-i, 
given  by  (1.9 )— ( 1.10)  with  n  =  m.  by  replacing  certain  node-weight  pairs.  This  can  be  carried 
out  by  successively  removing  an  node-weight  pair  by  Algorithm  2.1,  and  then  adding  a  new 
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node- weight  pair  u^inj;  one  stop  of  A lu,or it  Inn  '.i.l  i;i  i_s].  This  combination  of  algorithms  yields 
a  sliding  window  scheme. 

3      Downdating  and  the  FFT  algorithm 

When  the  nodes  Zk  are  equidistant  on  the  unit  circle  and  all  weights  wj.  are  unity,  then 
interpolation  and  least-squares  approximation  of  a  given  function  by  algebraic  or  trigonometric 
polynomials  can  be  carried  out  rapidly  by  the  FFT  algorithm.  This  section  describes  in  some 
detail  how  our  downdating  scheme  may  be  combined  with  the  FFT  algorithm  to  achieve  rapid 
schemes  for  interpolation  when  the  nodes  are  essentially  equidistant.  More  precisely,  we  will 
consider  the  case  when  the  set  of  nodes  {~a.}/."=i  '"  r'10  hitler  product  (1.1)  is  a  subset  of  the 
set  of  equidistant  nodes  {uxp{'2-jj/N )})St)\  where  k  :=  .V  -  m  is  a  small  positive  integer.  In 
our  operation  count  we  will  assume  that  n  is  independent  of  in.  The  weights  »•£  in  (  1.1)  are  all 
assumed  to  be  unity. 

Let  /'  denote  a  complex-valued  lunci  ion.  im-e  values  /(-/,-),  1  '_  />'  <  m.  are  explicitly 
known.  We  remark  that,  a  representation  ol  i  he  interpolating  polynomial  /Jm_i  £  I [m_ j  in 
Newton  form  can  be  computed  in  0(m2  j  arithmetic  operations.  The  Vandcrmonde  solver  by 
Bjorck  and  Percyra  [2]  can  be  used  to  determine  a  representation  of  }>,n-\  in  terms  of  the 
monomial  basis,  and  requires  also  O(m-)  arithmetic  operations.  Our  scheme  only  requires 
0(m  log  m)  arithmetic  operations  and  yields  a  representation  of  ;.>m_[  in  the  basis  of  Szego 
polynomials. 

Let  {-[.}£_[  denote  the  complement  of  the  set  {-(.Jl'Li  in  {expi  2~ij/.\T  ))^S^ .  and  let 
/{={.)  :=  0.  1  <  /.:  <  k.  Thus,  /is  defined  at  the  .Y  roots  of  unity  cxp(27r/(j  -  I  )/.Y).  I  <j  <  .V, 
and  we  can  compute  the  Fourier  coefficients  of  the  polynomial  /-».v-t  €  II.v-i  that  interpolates 
/  at  these  nodes  by  the  TFT  algorithm  in  0{\  log  -Y  )  =  0(m  log  m)  arithmetic  operations. 
Using  the  Schur  parameters  given  in  Example  1.1.  wo  apply  Algorithm  '_'.]  k  times  to  eliminate 
the  nodes  :'k,  1  <  k  <  k,  from  the  inner  product.  This  requires  O(m)  arithmetic  operations 
and  yields  the  Fourier  coefficients  of  the  desired  interpolating  polynomial  pm-\.  The  Fourier 
coefficients  of  the  least  squares  approximates  /;„_i  c  n„_[  with  n  <  m  are.  of  course,  a  subset 
of  the  Fourier  coellicieuts  ol  i>,,,-\- 

A  scheme  closely  related  to  the  one  outlined  above  can  be  used  to  compute  trigonometric 
polynomials  rapidly,  hot  :^  and  :'k  be  the  nodes  introduced  above,  and  define  Bk  :=  arg(r^.), 
!</,•<  in,  and  H'k  :=  arg(r[J.  !</.'<  n.   Also  assume  that  m  is  odd.  i.e..  in  =  2r  4-  1.   Lot 
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f{8)  be  a  real- valued  function  defined  at  the  nodes  #;.  and  let  f{9[)  :=  0.  1  <  A;  <  k.  We  wish 

to  compute  a  trigonometric  polynomial  i(9)  t  span{l.cos  9 cos  rtf.sin  9 sin  r9)  that 

minimizes  the  sum 

m 

£j/(^)-  ((9k))2  ■ 
fc=i 

This  can  be  accomplished  as  follows.  First  solve  the  minimization  problem 

min     J2]e2*i(k-l)r/Nf  f(l._  1}^  _£0ie2ri(*-l)(i-l)/Ar,2 

*,eC  fc=1  V  -A  /      J=i 

by  the  FFT  algorithm  in  0(m  log  m)  arithmetic  operations,  and  denote  the  optimal  coefficients 
by  aj.  Remove  the  nodes  9'k,  1  <  k  <  k,  by  applying  the  downdatinc;  scheme  k  times  to  the 
polynomial  /*2r+«(~)  :=  HjLu*  'V7'  where  A'  -  1  =  'J/'  +  k.  This  requires  O(m)  arithmetic 
operations  and  yields  the  polynomial  p?,-  €  ft?,-  'lie  desired  trigonometric  polynomial  is  then 
given  by  t{9)  :=  z~rp2r{z),  z  =  cxp(i'0). 

4      Computed  examples 

We  now  present  the  results  of  some  numerical  experiments  with  our  downdating  procedure. 
These  experiments  were  performed  in  FORTRAN'  on  a  SparcStation  SLC  at  Northern  Illinois 
University,  on  which  there  are  approximately  7  and  Hi  significant  decimal  digits  in  single- 
precision  and  double-precision  calculations,  respectively. 

The  first  set  of  experiments  compares  the  accuracy  of  the  downdating  procedure  with  that  of 
the  updating  procedure  IUQR  described  in  [8].  Wo  input  N  unimodular  nodes  {zj}jLi,  positive 
weights  {w^}^=l,  and  complex  function  values  {f{zj)}^Ll.  For  any  positive  integer  m  <  N,  let 
am  denote  the  vector  of  Fourier  coeificicnts  of  the  solution  pm-\  of  (1.6)  with  n  =  in.  and  let 
gm  denote  the  vector  of  Schur  parameters  determined  by  the  inner  product  (1.1)  and  recurrence 
relations  (1.2)-(1.3). 

We  first  compute  vectors  a.v  and  g/v  using  an  implementation  of  the  IUQR  algorithm  in 
single-precision  arithmetic  [8].  We  then  repeatedly  apply  a  single-precision  implementation  of 
Algorithm  2.1  to  compute  am  and  gm  for  decreasing  values  of  m.  The  kth  application  of  Al- 
gorithm 2.1  removes  the  node  rv-t+i  from  the  inner  product  to  compute  the  solution  of  (1.6) 
with  n  :=  in  :=  .V  -  /,-.  After  each  downdating  step,  we  calculate  the  relative  error  in  am,  i.e., 
||am  -  am||>/||am||2,  and  the  error  in  gm,  i.e.,  ||g,„  -  gm||2,  where  ||x||-2  denotes  the  Euclidean 
norm  of  x  6  Cm.  We  also  solve  each  problem  of  order  m  using  the  IUQR  algorithm  in  single- 
precision  arithmetic  and  compute  the  resulting  errors.   The  results  of  the  IUQR  algorithm  in 
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double-precision  arithmetic  arc  used  as  exact  answers  in  the  error  calculations.  The  following 
tables  display  the  resulting  errors  for  the  downdating  procedure  (DD)  and  the  updating  proce- 
dure (  I'D).  In  all  of  the  experiments,  each  function  value  /(-;)  lias  its  real  part  and  imaginary 
part  randomly  generated  according  to  a  uniform  distribution  in  [  —  5,5]. 

We  first  choose  the  nodes  to  be  the  N  :=  300  roots  of  unity  Zj  :=  exp(2-;(j  —  1)/N), 
1  <  j  <  Ar,  and  uniform  weights  wj  :=  1.  Table  1  shows  the  errors  for  problems  of  order 
in  :=  iV  —  k  for  various  values  of  A.1.  Table  1  also  shows  the  results  with  the  same  choice  of  nodes 
and  weights,  except  that  the  nodes  are  permuted  in  a  random  way.  This  permutation  changes 
the  nodes  that  arc  deleted  as  well  as  the  order  in  which  the  nodes  are  added  in  the  updating 
procedure.  It  should  be  noted  that  the  errors  in  the  downdating  procedure  can  be  expected  to 
increase  as  k  increases.  Table  2  shows  that  similar  results  are  obtained  with  the  same  set  of 
nodes  and  randomly  generated  weights  in  (0.10). 

Table  3  shows  the  results  obtained  with  uniform  weights  and  .V  :=  300  nodes  zj  :  = 
exp(~/(j  —  1)/.V).  1  <  j  <  .V.  in  [0.-).  both  in  their  original  order  and  in  a  random  order. 
Here  again,  the  downdating  procedure  perform-  well. 

Table  -I  shows  the  results  when  the  initial  300  nodes  are  randomly  selected  points  on  the 
unit  circle.  In  this  example,  the  error  in  the  downdating  procedure  displays  a  sudden  increase 
at  k  =  10.  In  this  step  the  error  in  the  computed  downdated  weight,  \iv  —  tt)|,  was  greater  than 
10"'.  Observe  that  the  error  incurred  at  this  downdating  step  propagated  to  the  subsequent 
downdating  steps,  but  the  errors  in  a..,  and  g,„  seem  to  grow  gradually  as  k  increases,  i.e.,  as 
m  :=  .V  —  k  decreases.  Other  experiments  with  random  nodes  produced  similar  results.  It 
should  be  noted  thai  in  our  experiments,  a  large  error  in  the  computed  weight  did  not  always 
coincide  with  a  large  jump  in  the  errors  in  am  and  g,„.  It  is  clear  that  more  work  needs  to  be 
done  in  order  to  understand  the  numerical  aspects  of  the  updating  and  downdating  problems. 

Our  final  experiment  tests  the  accuracy  and  speed  of  the  procedure  described  in  Section 
3  for  downdating  the  IT  FT.  The  .'V-point  FFT  is  used  to  obtain  the  Fourier  coefficients  of  the 
interpolating  polynomial  ]>,\-\.  A  randomly  selected  set  of  nodes  is  then  removed  from  the 
inner  product  using  Algorithm  2.1.  This  experiment  was  run  with  a  radix- 2  FFT  subroutine  in 
single  precision  arithmetic  with  .V  :—  1021  and  .V  :=  2018.  Table  5  shows  the  computed  error 
alter  /•;  downdating  steps  for  various  values  of  k.  As  above,  we  use  the  results  of  the  1UQR 
algorithm  in  double  precision  arithmetic  as  exact  answers  for  error  checking.  We  also  display 
the  amount  of  CPU  time  required  by  the  FFT  with  k  downdating  steps  and  the  time  required 
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Table  1:  300  nodes  with  equispaced  arguments  in  [0,2");  uniform  weights. 


nodes  in 

order  ol  incrc 

using  arguii 

icni . 

nodes  m  random  order. 

relative  error  in  nnt 

error 

in  'Xm 

relative  error  in  a,„ 

error 

"1   g-n 

k 

DD 

ID 

DD 

I'D 

DD 

UD 

DD 

ED 

(J 

0.122E-03 

0.122E-03 

0.12GE-03 

0.12(5  E-03 

0.2S4E-04 

0.284E-01 

0.310E-04 

U.340E-01 

10 

0.1  LlE-03 

0.121E-03 

0.295E-03 

0.1  10  E-03 

0.335E-04 

0.290E-01 

0.545E-04 

0.340E-04 

20 

0.M1E-03 

0.115E-03 

0.405E-03 

0.H8E-03 

0.473E-04 

0.314E-04 

0.211  E-03 

0.361E-04 

30 

0.115E-03 

0.100E-03 

0.322E-03 

0.154  E-03 

0.550E-04 

0.272E-04 

0.272E-03 

0.331E-04 

-10 

0.112E-03 

0.967E-0-1 

0.340E-03 

0.1 55  E-03 

0.943E-04 

0.257E-04 

0.359E-03 

0.374E-04 

50 

0.112E-03 

0.940E-04 

0.329E-03 

0.164E-03 

0964E-04 

0.259E-04 

0.325E-03 

0.437E-04 

70 

0.124E-03 

0.102E-03 

0.578E-03 

0.171E-03 

0.174E-03 

0.258E-04 

0.424E-03 

0.422E-04 

90 

0.132E-03 

0.1 05  E-03 

0.725E-03 

0.173E-03 

0.254E-03 

0.235E-04 

0.677E-03 

0.422E-04 

110 

0.154  E-03 

0.107E-03 

0.305E-03 

0.104E-03 

0.375  E-03 

0.238E-04 

0.753E-03 

0.383E-04 

130 

0.151E-03 

0.117E-03 

0.2S9E-03 

0.160E-03 

0.465E-03 

0.226E-04 

0.109E-02 

0.359E-04 

150 

0.1G6E-03 

0.127E-03 

0.321  E-03 

0 .1 50  E-03 

0. GOO  E-03 

0.18GE-04 

0.119E-02 

0.463E-04 

by  the  single-precision  ll'QR  algorithm  on  the  problem  of  order  m  =  X  —  k.  It  is  interesting  to 
note  that  the  downdating  procedure  produces  substantially  more  accurate  answers  faster  than 
the  IUQR  algorithm  even  lor  moderately  sized  values  of  k. 

5      Conclusion 

New  algorithms  are  presented  for  discrete  least  squares  approximation  by  algebraic  poly- 
nomials on  the  unit  circle  and  by  trigonometric  polynomials.  The  algorithms  downdate  the 
approximant  as  well  as  recurrence  coefficients  for  Szcgo  polynomials  when  a  node-abscissa  pair 
is  removed  from  the  inner  product.  The  schemes  arc  based  on  the  QR  algorithm  for  unitary 
Ilcssenberg  matrices.  A  particularly  attractive  application  is  to  combine  our  algorithm  for 
trigonometric  polynomials  with  the  fast  Fourier  transform  algorithm.  This  yields  a  fast  scheme 
for  computing  trigonometric  least  squares  approximants  when  the  nodes  of  the  inner  product 
form  a  subset  of  equidistant  nodes  on  [0,2"). 
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Tabic  2:  300  nodes  with  equispaced  arguments  in  [0.2rr);  random  weights  in  (0.101 


nodes  in  ( 

niier  ol  incrc 

asing  iirgum 

int. 

nodes  in  random  order. 

relative  error  in  am 

error 

"I   Km. 

relative  error  in  am 

error 

"1   gm 

k 

DD 

ID 

DD 

I'D 

DD 

UD 

DD 

UD 

0 

0.207E-03 

0.207E-U3 

0.2-1-1  E-03 

0.2-HE-03 

0.511E-04 

0.511E-0-1 

0.572E-04 

0.572E-04 

10 

0.206E-03 

0.183E-03 

033")  E-03 

0.184  E-03 

0.563  E-04 

0.537E-04 

0.844E-04 

0.744E-0I 

20 

0.1 99  E-03 

0  171  E-03 

0.417E-03 

0.1 92  E-03 

0.029E-0-1 

0.512E-04 

0.140E-03 

0.722E-04 

:iU 

0.1S7E-03 

0.173E-03 

0.302E-03 

0.198E-03 

0.7S5E-04 

0.518E-04 

0  252E-03 

0.750E-04 

■10 

0.188E-03 

0.1 79  E-03 

0.340E-03 

0.203  E-03 

0.873E-04 

0.509E-04 

0.303E-03 

0.739E-04 

50 

0.181  E-03 

0.1G7E-03 

0.331E-03 

0.209  E-03 

0.30SE-03 

0.487E-04 

0.414E-03 

0.714E-01 

70 

0.206  E-03 

0.109  E-03 

0.604  E-03 

0  217E-03 

0.406E-03 

0.479  E-04 

0.024  E-03 

0.038E-04 

90 

0.192E-03 

0  11.")  E-03 

0.7G5E-03 

0.222  E-03 

0.572E-03 

0442E-01 

0.17SE-02 

0.088E-01 

110 

0.19SE-03 

0.155  E-03 

0.389  E-03 

0.220E-03 

0.S74E-03 

0.333E-04 

0.184E-02 

0.533E-01 

130 

0.201  E-03 

0.155E-03 

0.3-1 1  E-03 

0.213E-03 

0894E-03 

0.286E-04 

0.139E-02 

0.401E-04 

150 

0.201  E-03 

0.160E-03 

0.379E-03 

0.1 95  E-03 

0.1 22  E-02 

0.256E-04 

0.153E-02 

0.2S5E-04 

Table  3:  300  nodes  with  equispaced  arguments  in  [0.x);  uniform  weights. 


nodes  in  order  of  increasing  argument. 

nodes  m  random  order. 

relative  error  m  a,„ 

error 

i"  g». 

relative  error  in  am 

error 

">  g»» 

k 

DD 

UD 

DD 

UD 

DD 

UD 

DD 

UD 

0 

0.398E-03 

0.398  E-03 

0.7  12  E-03 

0.742E-03 

0.125E-03 

0.1 25  E-03 

0.110E-03 

0.116E-03 

10 

0.455E-03 

0.398E-03 

0.188E-02 

0.743  E-03 

0.125E-02 

0.101E-03 

0.106E-02 

0.124E-03 

20 

0.437E-03 

03S9E-03 

0.1S7E-02 

0.720E-03 

0.1 44E-02 

0.112E-03 

0.146E-02 

0.122E-03 

30 

0.452E-03 

0.390E-03 

0.22SE-02 

0  71 4  E-03 

0170E-02 

0.110E-03 

0.163E-02 

0.131E-03 

40 

0.479E-03 

0.391  E-03 

0.205E-02 

0.694  E-03 

0.195E-02 

0.947E-04 

0.1 67  E-02 

0.113E-03 

50 

0.508E-03 

0.391  E-03 

0.271  E-02 

0.680E-03 

0  215E-02 

0.886E-0  1 

0.1 65  E-02 

0.1 09  E-03 

70 

0.57GE-03 

0.131E-03 

0.278E-02 

0.045E-03 

0.251  E-02 

0.796E-04 

0.210E-02 

0.1 09  E-03 

90 

0.533E-03 

0.150E-03 

0.1 08  E-02 

0.605E-03 

0  367 E-02 

0.849E-04 

0.302  E-02 

0.1 06  E-03 

110 

0.558E-03 

0.435E-03 

0  193  E-02 

0.552  E-03 

0.158E-02 

0.795E-01 

0.367E-02 

0.102E-03 

130 

0.650  E-03 

ii  l23E-():i 

0  2 10  E-02 

0.505E-03 

0.528  E-02 

0.792E-OI 

0.539E-02 

O.S40E-0I 

1 50 

0.110E-02 

1)  I93E-03 

il  350E-02 

n   1(56  E-03 

0.671  E-02 

0.548E-04 

0.528E-02 

0.895E-0  1 
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Tabic  -1:  300  nodes  with  random  arguments  in  [0,2~);  uniform  weights. 


I'lative  e 

■ror  in  ;»,„ 

error 

">  gm 

k 

DD 

I'D 

DD 

UD 

i) 

n  i  in:-03 

0  1 1LE-03 

0.025  E-03 

0.925E-03 

9 

0.50SE-03 

0.  IOOE-03 

O.OME-03 

0.S38E-03 

10 

U.2-18E-0I 

0.1G0E-03 

0.31SE-01 

0.S38E-03 

20 

0.212E-01 

0.129E-03 

0.373E-01 

0.789E-03 

30 

0.220E-01 

0  KME-03 

0  32SE-01 

0.735E-03 

•10 

0.1ME-01 

0.117E-03 

0.131E-01 

0.826E-03 

:>() 

0.18GE-01 

0.124  E-03 

0.240E-01 

0.654E-03 

70 

0.106E-01 

0.137E-03 

U.2G1E-01 

0.658E-03 

00 

0.210E-01 

0.122E-03 

0.275E-01 

0.483E-03 

110 

0.3  17E-01 

0.137E-03 

0.820E-01 

0.100E-02 

130 

0.502E-01 

0.120E-03 

0.852E-01 

0.508E-03 

150 

0.513E-01 

0.7S8E-04 

0.865E-01 

0.3G0E-03 

Lfj 


Table  5:  Downdating  the  FFT 


.V  = 

in'.'l 

relative  error  in  a,„ 

error 

in  g,n 

CPU  seconds 

in 

it 

DD 

I'D 

DD 

UD 

DD 

UD 

1013 

11 

0.4G1E-05 

0.213E-01 

0.21SE-05 

0.204E-01 

0.303E+01 

0.148E+03 

993 

31 

U.502E-05 

0.215E-01 

0.623E-05 

0.195E-01 

0.814E+01 

0.142E+03 

973 

51 

0.603E-05 

0.217E-01 

0.142E-04 

0.193E-01 

0.132E+02 

0.136E+03 

923 

101 

0.1)00  E-05 

0.193E-01 

0.257E-04 

0.244E-01 

0.253E+02 

0.122E+03 

913 

111 

0.944E-05 

0.191E-01 

0.2G5E-04 

0.242E-01 

0.276E+02 

0.119E+03 

S93 

131 

O.lllE-0'1 

0.17SE-01 

0.312E-04 

0.219E-01 

0.323E+02 

0.114E+03 

873 

151 

0.129E-04 

0.174E-01 

0.373E-04 

0.220E-01 

0.368E+02 

0.109E+03 

853 

171 

0.143E-04 

0.172E-01 

0.462E-04 

0.221E-01 

0.412E+02 

0.104E+03 

833 

191 

0.172E-04 

0.174E-01 

0.566E-04 

0.220E-01 

0.456E+02 

0.990E+02 

813 

211 

0.223E-04 

0.171E-01 

0.721E-04 

0.241E-01 

0.497E+02 

0.941E+02 

713 

311 

ii  177E-0! 

D.133E-01 

0.131E-03 

0.169E-01 

0.692E+02 

0.719E+02 

G13 

■111 

0.167E-03 

0.104E-01 

0.458E-03 

0275E-01 

0.862E+02 

0.527E+02 

513 

"ill 

0.2G7E-03 

U  7G9E-02 

0.527E-03 

0.174E-01 

0.101E+03 

0.365E+02 

■113 

Gil 

0.28GE-03 

0.417E-02 

0.487E-03 

0.852E-02 

0.112E+03 

0.233E+02 

313 

711 

0  337E-03 

0.322E-02 

0.55-1  E-03 

0.781E-02 

0.122E+03 

0.132E+02 

.V  = 

2048 

relative  error  in  a,„ 

error 

in  gm 

CPU  seconds 

m 

k 

DD 

UD 

DD 

UD 

DD 

UD 

2037 

11 

0.598E-05 

0.105E+00 

0.259  E-05 

0.184E  +  00 

0.391E+01 

0.390E+03 

1997 

51 

0.740E-05 

0.110E+00 

0.151E-04 

0.235E+00 

0.172E+02 

0.374E+03 

1947 

101 

0.100E-04 

0.102E+00 

'|  293E-04 

0.164E+00 

0.336E+02 

0.355E+03 

1897 

151 

0.125E-04 

0.980E-01 

0.426E-04 

0.126E+00 

0.495E+02 

0.337E+03 

1847 

201 

0.159E-04 

0.9ISE-01 

0.G42E-04 

0.131E+00 

0.649E+02 

0.319E+03 

1797 

251 

0.237E-04 

0.881E-01 

0.812E-04 

0.136E+00 

0.799E+02 

0.301E+03 

1747 

301 

i)  289E-0I 

i»  8G5E-01 

0.1 10  E-03 

0.157E+00 

0.94GE+02 

0.2S5E+03 

1097 

35 1 

0.358E-0  1 

0  834E-01 

0.1 15  E-03 

0.137E+00 

0.109E+03 

0.26SE+03 

1647 

•101 

0.514E-04 

0  809E-0I 

0  170  E-03 

0.125E+00 

0.123E+03 

0.252E+03 

1597 

151 

0.749E-0I 

0.748E-01 

0.2I2E-03 

0.133E+00 

0.136E+03 

0.237E+03 

1547 

501 

0  779E-04 

0.729E-01 

0  2  13  E-03 

0.132E+00 

0.149E+03 

0.222E+03 

1447 

(301 

0  953E-0  1 

0.G95E-01 

0  313E-03 

0.149E+00 

0.174E+03 

0.194E+03 

1347 

701 

0.242E-03 

0.G28E-01 

0  39GE-03 

0.114E+00 

0.197E+03 

0.167E+03 

1147 

901 

0.517E-03 

0143E-01 

0.109E-02 

0.116E+00 

0.238E+03 

0.121E+03 

1047 

1001 

0.731E-03 

0.381E-01 

0.113E-02 

0.640E-01 

0.25GE+03 

0.100E+03 

947 

1101 

0.114E-02 

0.285E-01 

0.183E-02 

0.523E-01 

0.273E+03 

0  815E+02 
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